library(MCMCpack)
library(foreign)
library(car)

save <- TRUE
fig.size <- 12

setwd("~/Google Drive/Broockman/Completed/Extremity vs Consistency/ReplicationData")

#Get mass data
mass.data <- read.dta("MainSurveyApril2013/ssi-survey-apr-2013.dta", convert.underscore = TRUE)
names(mass.data)

#Get legislator data
legislator.data <- read.dta("StateLegislatorSurvey/state-legislator-survey.dta", convert.underscore = TRUE)
names(legislator.data)



#######MULTIS#########

#####POLICY LEVEL
#Mass Public
mass.multis <- mass.data[,22:33]
mass.extremity <- abs(mass.multis - 4)
mass.extremity <- rowMeans(mass.extremity, na.rm = TRUE)

#Elites
elite.multis <- legislator.data[,21:32]
elite.extremity <- abs(elite.multis - 4)
elite.extremity <- rowMeans(elite.extremity, na.rm = TRUE)

if(save) pdf("MiscFigures/DATA_mainfig_typical_extremity_of_responses.pdf", width = fig.size, height = fig.size, pointsize = 16)
plot(NA, main = "Typical Extremity of Subjects' Responses",
     ylab = "Density", xlab = "Typical Extremity of Subjects' Responses",
     xlim = c(0.3,2.5), ylim = c(0,1.3), yaxt='n')
lines(density(mass.extremity, na.rm = TRUE, weights = mass.data$weights), lwd = 2, lty = 1)
lines(density(elite.extremity, na.rm = TRUE), lwd = 2, lty = 4)
text(2,.8,labels="Mass Public", cex = 1.3)
text(.75,1,labels="Legislators", cex = 1.3)
if(save) dev.off()

###RESPONDENT LEVEL SUMMARIES
mass.issue.specific.averages.across.rs <- abs(rowMeans(mass.multis, na.rm = TRUE) - 4)
elite.issue.specific.averages.across.rs <- abs(rowMeans(elite.multis, na.rm = TRUE) - 4)

#Respondent Level
if(save) pdf("MiscFigures/DATA_mainfig_respondents_average_position.pdf", width = fig.size, height = fig.size, pointsize = 16)
plot(NA, main = "Extremity of Subjects' 'Average Response' Across Issues", ylab = "Density", xlab = "Extremity of Subjects' 'Average Response'", xlim = c(0,2.5), ylim = c(0,1))
lines(density(mass.issue.specific.averages.across.rs, na.rm = TRUE, weights = mass.data$weights), lwd = 2, lty = 1)
lines(density(elite.issue.specific.averages.across.rs, na.rm = TRUE), lwd = 2, lty = 4)
text(0.3,.9,labels="Mass Public", cex = 1.3)
text(1.3,.8,labels="Legislators", cex = 1.3)
if(save) dev.off()


####IRT
#Isolate IRT questions
binary.iqs.mass <- mass.data[,c(1:9,11:21)]
binary.iqs.legislators <- legislator.data[,1:20]


#Remove legislators who didn't answer
count.nas <- rep(0,times=nrow(legislator.data))
for(i in 1:length(binary.iqs.legislators)){
	count.nas <- count.nas + is.na(legislator.data[,i])
}
table(count.nas)
binary.iqs.legislators <- binary.iqs.legislators[which(count.nas < 5),]


#Combine datasets 
cbind(names(binary.iqs.mass), names(binary.iqs.legislators)) #check that they line up
all.iqs <- rbind(binary.iqs.mass, binary.iqs.legislators)
is.legislator <- c(rep(0,times=nrow(binary.iqs.mass)), rep(1,times=nrow(binary.iqs.legislators)))


#Build ideology IRT model from IQ responses
typical.ideology.scale <- colMeans(MCMCirt1d(all.iqs))

#Look at question response distributions
if(FALSE){
	typical.ideology.scale.items <- colMeans(MCMCirt1d(binary.iqs, store.ability = FALSE, store.item = TRUE))
	typical.ideology.scale.items
	colMeans(binary.iqs, na.rm = TRUE)
}

if(save) pdf("MiscFigures/DATA_all_together_now.pdf", width = fig.size, height = fig.size, pointsize = 16)
plot(NA, main = "Distribution of Cross-Policy 'Ideal Points'", ylab = "Density", xlab = "", xlim = c(-3.5,3.5), ylim = c(0,.7))
lines(density(typical.ideology.scale[which(is.legislator==0)], na.rm = TRUE), lwd = 2, lty = 1)
lines(density(typical.ideology.scale[which(is.legislator==1)], na.rm = TRUE), lwd = 2, lty = 4)
text(1.4,.5,labels="Mass Public", cex = 1.3)
text(-2.7,.25,labels="Legislators", cex = 1.3)
if(save) dev.off()